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1.  Introduction. 

In  computing  the  variance  of  a  sample  of  N  data  points  {x^},  the  fundamental 
calculation  consists  of  computing  the  sum  of  squares  of  deviations  from  the  mean. 
This  quantity,  which  for  brevity  will  be  referred  to  as  “the  sum  of  squares"  or 
simply  as  “5” ,  is  defined  as 


where 


N 


s  ~~  x)2> 
1=1 


(1.1a) 


(1.16) 


This  computation  can  be  easily  performed  directly  by  the  two-pass  algorithm 
(1.1)  provided  that  (a)  N  is  small  compared  to  the  amount  of  core  memory  avail¬ 
able  and  (b)  the  variance  is  not  too  small  relative  to  the  norm  of  the  data, 
||x||2  =  x^)1/2.  If  either  of  these  conditions  is  violated,  however,  the  situa¬ 

tion  changes.  If  N  is  large,  this  approach  to  computing  S  may  be  unsatisfactory 
since  it  requires  passing  through  the  data  twice:  once  to  compute  x  and  then  again 
to  compute  S.  This  problem  is  sometimes  avoided  by  use  of  the  following  textbook 
algorithm,  so  called  because,  unfortunately,  it  is  often  suggested  in  statistical 
textbooks: 


p-m 


(2) 


This  rearrangement  allows  S  to  be  computed  with  only  one  pass  through  the  data, 
but  the  computation  may  be  numerically  unstable  and  should  almost  never  be 
used  in  practice.  This  instability  is  particularly  troublesome  when  S  is  very  small 
compared  to  ||x||2,  in  which  case  even  the  two-pass  algorithm  can  be  unstable. 

In  discussing  the  stability  of  a  numerical  scheme  for  computing  S,  a  useful 
concept  is  that  of  the  condition  number  k  of  the  data.  This  quantity  was  first 
introduced  by  Chan  and  Lcwis[2]  who  give  a  thorough  discussion.  Briefly,  k  is  a 
measure  of  the  sensitivity  of  S  to  changes  in  the  data.  The  quantity  ku  is  an  upper 
bound  for  the  relative  perturbation  which  would  occur  in  the  exactly  computed 
sum  of  squares  if  the  input  data  contained  relative  errors  of  size  u.  If  the  true  sum 
of  squares  is  S,  then  k  is  given  by 


K  = 


IkHa 

Vs ' 


(1.3) 
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It  is  easy  to  see  that  *  >  1  and  that  in  general  k  grows  as  the  variance  decreases. 

An  error  analysis  of  the  textbook  algorithm[2]  shows  that  the  relative  error 
in  5  can  be  bounded  by  something  on  the  order  of 

3  Nk2u, 

where  u  is  the  machine  roundoff  unit  (see  section  6).  This  algorithm  is  therefore 
seldom  useful,  as  confirmed  by  the  experimental  results  of  Table  1. 

The  error  analysis  of  the  two-pass  algorithm  found  in  section  6  shows  that 
the  relative  error  in  the  sum  of  squares  computed  using  that  algorithm  cun  be 
bounded  by 

Nu  -f-  N2k2u 2. 

The  second  term  in  this  bound  has  traditionally  been  ignored  in  error  analyses  of 
the  two-pass  algorithm  as  being  of  second  order.  But  in  the  ease  we  are  interested 
in  here,  when  N  and  k  arc  both  large,  this  term  can  easily  dominate.  Table  2 
shows  this  happening  in  practice. 

During  the  preparation  of  this  manuscript,  a  simple  modification  of  the  two- 
pass  algorithm  was  found  by  Professor  Ake  Bjorck  which  reduces  this  bound. 
Based  on  the  error  analysis  of  section  6  for  the  standard  two-pass  algorithm,  he 
suggested  computing'5  by 

N  fN  \  2 

s -  w(£(l' ~  i)J  ■  <*• 

In  exact  arithmetic  the  second  term  is  zero,  but  computationally  it  is  a  good 
approximation  to  the  error  in  the  two-pass  algorithm.  Note  that  (1.4)  can  also 
be  viewed  as  the  textbook  algorithm  applied  to  the  data  {(x, —  x)}.  The  error 
analysis  of  section  6  shows  that  the  relative  error  in  S  computed  by  (1.4)  can  be 
bounded  by 

Nu  +  4N2«u2. 

This  modification  adds  only  N  additions  and  2  multiplications  to  the  cost  of  the 
two-pass  algorithm  (already  3 N  —  2  additions  and  /V-f-  1  multiplications)  and  can 
be  very  useful  when  the  data  is  poorly  conditioned.  See  table  3  for  some  numerical 
results. 

Of  course  formula  (1.4)  is  still  a  two-pass  algorithm.  For  large  N  it  may 
be  desirable  to  compute  5"  with  only  one  pass  through  the  data.  A  number  of 
papers  have  appeared  recently  on  "updating"  algorithms  for  computing  S.  These 
are  algorithms  which  are  based  on  formulae  for  adding  one  new  data  point  to  a 


sample  and  computing  the  value  of  S  for  the  combined  sample  by  updating  the 
(presumably  known)  value  of  S  for  the  original  sample.  By  starting  with  a  sample 
of  size  1  and  applying  this  formula  repeatedly,  we  get  a  onc-pnss  algorithm  for 
computing  S  for  a  sample  of  arbitrary  size.  Youngs  and  Cramcr[G)  have  inves¬ 
tigated  several  such  algorithms  and  have  found  the  following  algorithm  to  be  the 
best: 


S  :=  0 
T x\ 

for  j  :=  2, 3, . . N  do 
T:=T  +  x. 


S:=S  + 


1 

Hi  - 1) 


Uxj 


(1.5) 


This  is  based  on  the  updating  formula 


Sij  —  + 


1 

;(;  —  i) 


(IB) 


where  Sij  stands  for  the  sum  of  squares  for  the  data  points  through  x3  and  TltJ 
is  the  sum  of  through  Xj.  This  notation  will  be  used  throughout. 

One  imporant  characteristic  of  this  updating  formula  is  that  S\j  is  formed 
from  Sij—i  by  adding  to  it  a  nonnegative  quantitiy.  In  the  textbook  Algorithm 
(1.2),  on  the  other  hand,  S  is  formed  by  a  subtraction  which  can  lead  to  gross 
cancellation  and  even  to  negative  values  of  S  being  computed. 

In  practice  the  method  (1.5)  generally  performs  on  a  level  comparable  to  the 
two-pass  algorithm  (1.1).  Chan  and  Lewis[2]  present  detailed  error  annly:.  >  of 
some  similar  updating  methods. 

In  the  next  section  we  present  a  generalization  of  the  updating  formula  (l.C) 
for  combining  two  samples  of  arbitrary  size.  Then  in  section  3  we  describe  a 
pairwise  algorithm  for  computing  S  which  is  essentially  still  a  one-pass  algorithm 
but  which  numerically  is  often  more  stable  than  the  standard  two-pass  algorithm. 


2.  A  General  Updating  Formula. 

The  method  of  Youngs  and  Cramer  depends  on  an  updating  formula  which 
allows  one  to  computed  forj-j- 1  points  when  given  the  value  of  S  for  j  |K>iius  and 
one  new  point.  In  other  words,  we  can  combine  a  sample  of  size  j  with  a  sample 
of  size  1  and  determine  the  value  of  5  for  the  combined  sample. 

This  formula  can  be  easily  generalized  to  allow  us  to  combine  two  samples  of 
arbitrary  size.  Suppose  we  have  two  samples  {x,-}£Lj,  and  v.c  know 


3 


m 


m-(-n 


^m+l,m+n  “  ^  i  **» 

i»=l  i=m+l 


»=1 


m  .  m-(-n 

“  ^J(*»  ^1  ,tn)^»  “Sm-t-Um+n  =  V  ,  (®»  Tmq-](m-|  fl)  • 

m  .  ,  n 

i=rn-(-] 


Then,  if  we  combine  all  of  the  data  into  a  sample  of  size  m  -f-  n,  it  can  be  shown 

that 


Tj ,rn-f  n  “  Tiifn  -J-  Tj^-f-n 

5l,m+n  *=  ,m  “t"  'S’m+l.m-f-n 

+  n(m  +  n)(mT''"  ~  Tm+I'’"+") 


(2.1a) 

(2.16) 


If  we  rewrite  the  latter  formula  as 

Si,wi-f  n  —  Sj  m  S'm+l.m+n  *f* 


m 


n(m  -f 


then  wc  see  that  for  m  =  1,  n  *=  j  —  1,  this  reduces  to  the  formula  of  Youngs 
and  Cramer,  since  S  *=*  0  for  any  single  data  point.  The  form  (2.1b)  is  more  stable 
numerically,  however. 

Regardless  of  what  method  is  used  to  compute  S%  the  formulae  (2.1)  may  be 
useful  in  their  own  right  whenever  two  samples  must  be  combined.  One  possible 
application  is  to  parallel  processing.  If  one  has  two  or  more  prorr:  :ors  available,  the 
sample  can  be  split  up  into  smaller  subsamplcs,  and  the  sum  of  squares  computed 
for  each  subsamplc  independently  using  any  algorithm  desired.  The  sum  of  squares 
for  the  original  sample  can  then  be  calculated  using  the  updating  formal; 

However,  even  in  the  case  of  a  single  processor,  it  is  very  desirable  to  compute 
S  using  (2.1).  The  method  (1.5)  may  be  generalized  to  compute  S  by  processing 
the  data  in  groups  of  »n  elements:  compute  the  sum  of  squares  for  each  group  using 
the  two- pass  algorithm  and  then  update  the  global  S  accordingly.  Traditional 
updating  algorithms  such  as  that  of  Youngs  and  Cramer  have  used  w  1. 

We  have  found,  however,  that  the  stability  of  the  algorithm  is  increased  by- 
taking  m  >  1.  One  can  easily  sec  that  the  total  number  of  arithmetic  operations 
performed  on  the  data  is  minimized  by  taking  m  ***  \/N.  Wc  might  expect  that 
this  choice  of  m  will  minimize  the  resulting  error.  Although  >vc  do  not  have  a 
satisfactory  error  analysis  of  this  algorithm,  the  experimental  results  of  table  <1  do 
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tend  to  confirm  this  prediction.  Strictly  speaking,  with  m  >  1  this  is  no  longer  a 
one-pass  algorithm,  but  we  sec  that  only  m  data  values  at  a  time  need  to  be  kept 
in  core,  and  m  can  be  as  small  as  necessary. 

3.  The  Pairwise  Algorithm. 

Table  4  shows  that  choosing  m  >  1  not  only  gives  more  accuracy  than  using 
m  =  1,  but  can  actually  give  significantly  more  accuracy  than  the  tv. o- pass  algo¬ 
rithm.  This  suggests  that  when  computing  the  sum  of  squares  for  the  subsnuiplcs 
of  size  m,  we  should  not  use  the  two-pass  algorithm  when  £  is  small.  Hal  her,  \vc 
should  split  the  subsamplc  into  yet  smaller  groups.  Taking  this  idea  to  the  limit 
yields  a  pairwise  algorithm  analogous  to  the  well-known  pairv, ;so  algorithm  for 
computing  the  6um  of  N  numbers.  Let  5,(J  stand  for  the  sum  of  squares  of  elements 
Xi  through  Xj,  and  let  m  —  j_A//2j,  the  largest  integer  not  exceeding  N /2.  Then 
the  method  consists  of  computing  £i,n  by  first  computing  £]>m  and  £,„  \  \,w  and 
then  combining  these  by  means  of  (2.1).  Each  of  these  latter  quantities  has  been 
computed  by  a  similar  combination  of  still  smaller  subsamplcs. 

The  algorithm  can  be  implemented  as  just  described,  but  for  reasons  \.!.ich 
we  will  explain  shortly  it  is  actually  best  to  perform  the  pairwise  algorithm  in 
a  somewhat  modified  manner.  Consider  the  following  example  with  /V  =  13. 
Schematically,  we  compute  from  left  to  right  in  the  tableau  (3.1).  The  intermediate 
Ti,i  are  also  computed  in  a  similar  tableau  for  use  in  updating  the  Si  y  The  final 
value  TytN  will  be  the  sum  of  all  the  data  points  as  computed  by  the  pairwise 
summation  algorithm.  In  practice  we  can  compute  from  top  to  bottom  in  these 
tableaux  requiring  only  one  pass  through  the  data  and  using  only  0( log?  N)  si  orage 
locations  for  intermediate  results.  We  require  one  such  location  for  each  column 
in  each  tableau.  The  computation  for  the  tableau  (3.1)  would  proceed  as  follows: 

(a)  Compute  £1,2  and  store  in  Sjl], 

(b)  Compute  £3,4,  combine  with  S(l]  to  get  £,4,  and  store  this  in  S[2). 

(c)  Compute  5s(8  and  store  in  S[l]. 

(d)  Compute  S^g,  combine  with  S[l]  to  get  £5,8  and  then  combine  this 

with  S[2]  to  get  £i(8,  which  is  then  stored  in  S[3). 

(e)  Compute  <Sg,io  and  store  in  S[l]. 

(f)  Compute  £11,12,  combine  with  S[l]  to  get  S9  12,  which  is  then  stored 

in  S[2]. 

(g)  Clean-up  (necessary  when  N  is  not  a  power  of  2): 

Combine  X13  with  S{2]  to  get  69,13.  Combine  this  with  S[3]  to  get  £1,13. 
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Alternatively,  we  can  use  a  stack  structure  for  the  temporary  locations  ns  is  done 
in  the  sample  FORTRAN  routine  given  in  section  9. 

The  final  step  (g)  of  our  algorithm  requires  the  combination  of  samples  of  quite 
disparate  sizes.  Such  a  calculation  would  be  avoided  if  we  adopted  the  algorithm 
as  described  in  the  first  paragraph  of  this  section.  For  the  pairwise  summation 
algorithm,  Tsao[4]  points  out  that  the  corresponding  method  gives  a  decreased 
"average  error  complcxitj'"  and  presents  an  implementation  based  on  the  binary 
expansion  of  N.  That  strategy  could  be  adopted  in  the  present  context  as  well, 
but  its  usefulness  here  is  questionable.  We  feel  that  the  small  increase  in  accuracy 
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which  might  result  would  be  more  than  offset  by  the  increased  work  which  wc 
would  thus  incur.  For  the  updating  formula  (2.1b),  it  is  desirable  to  have  n  —  m 
whenever  possible,  since  that  formula  then  becomes  simply 


&1,2 m  —  *Si,m  ■5rn+J,2m  "h  T'm-\- i ,2 tn)  • 


In  this  respect,  the  tableau  (3.1)  gives  the  preferable  computational  scheme.  In 
fact,  the  amount  of  work  required  to  perform  the  pairwise  algorithm  ns  described 
here  is  not  significantly  more  than  that  required  for  the  two-pass  algorithm.  An 
operation  count  shows  that  roughly  2 N  additions  and  5N/2  multiplications  arc 
required,  as  opposed  to  3 N  additions  and  N  multiplications  for  the  tv.o-pass  al¬ 
gorithm.  In  addition,  some  bookkeeping  operations  are  required  to  manage  the 
pairwise  algorithm. 

Although  we  are  not  able  to  provide  any  error  bounds  proving  the  superiority 
of  the  pairwise  algorithm,  our  experimental  results  have  been  quite  satisfactory. 
Some  of  these  results  are  shown  in  table  5  of  section  8. 


5.  Extensions. 

Often  one  wants  to  compute  a  weighted  sum  of  squares  of  deviations  from 
the  mean, 


The  updating  formulae  (2.1)  still  hold  with  only  a  few  minor  modifications.  Let 
Wj'k  —  wi-  Then  (2.1)  is  replaced  by 

■  =  Tiifn  -J- 

cAW)  __  c (W)  j _ _ 

l,m+n  l,m  '  Wm+ltin+n(Wltin  +  Wm+1,m.,.M)  (5.2) 


x{  w,.„  r,'m 


Another  quantity  which  is  often  of  interest  is  the  covariance  of  two  samples 
{xi}  and  {y%}.  For  this  it  is  necessary  to  compute 


Ci,N  —  —  I)(y<  —  y). 
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If  we  let  Tfy  —  T’jfi  —  £*_yV»i  then  the  updating  formula  for  C  is 


^i,m+n  —  Cl,m  ~f"  C>n-)-l(m+n  ~}“ 


n(m  +  n) 


x  ( UtK*)  _ y(*)  \(  H.T(y)  _ r(y)  ^ 

■*  m-fl,m+n  M  my  l,m  1  m-j-l,m+nj- 


6.  Error  analysis  of  the  two-pass  algorithms. 

We  assume  throughout  our  error  analyses  that  wc  are  dealing  with  a  machine 
with  a  guard  digit  and  relative  precision  u.  On  a  base  /?  machine  with  a  t  digit 
mantissa  and  proper  unbiased  rounding,  u  =  *. 

Homan  letters  with  tildas  over  them  will  be  used  to  denote  quantities  ac¬ 
tually  computed  numerically.  The  same  letter  without  a  tilda  will  indicate  the 
corresponding  exact  quantity. 

In  this  section  we  present  error  analyses  for  both  the  standard  two-pas.-,  algo¬ 
rithm  (1.1)  and  Bjorck’s  modification  (1.4).  Let 

N 

Si  —  £](*,  —  i)2, 

»=i 

5  =  Si  -  Si. 

The  standard  two-pass  algorithm  is  5j.  We  first  compute  a  value  x  for  the 
mean  of  {x,}.  If  this  is  computed  in  the  standard  manner  wc  have 

,  N 

i  =  ±  V! *i(i  +  6),  with  161  <Nu  +  o(»J), 

N  ,t1  (0.1) 

The  computed  value  5j  is  then  given  by 
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Si 


—  Vi*  —  i)2(1  +  'h)>  I>J.|  <  (N  4-  2)u  +  0(u2) 

=  2  ((x*  —  *)  +  (*  —  *))  C1  +  %) 

=  Xj  ((*»  —  f)2  +  2(I»  ~  *)(*  ~  i)  +  (i  —  i)2)(l  4~  Vi) 

=  5  +  X](x*  ~  i)2,?‘  ~  (2  *•&)  2(J‘  ”  J)(1  +  7h) 


+(*2>*)V+x)4 


(0.2) 


The  0(u2)  terms  in  the  bounds  for  |r/,j  and  |£,|  turn  out  to  be  unimportant  in  the 
present  error  analysis  and  will  be  dropped  below.  Note  that  )T)(x,  —  x)  =  0  and 
that  the  following  inequalities  hold: 


|X>,  -  *)  Vj  <  Sltolloo  <  S(N  +  2 )u, 

|2>&|  <  W*hMh  <  N1/2||*yieiico  <  N3/2||x||2u, 

ID*.-  -  i)Vi  <  S1/2||r?||2  <  S^N1'2]^  <  S1/2A^2( A'  +  2)u. 


So  from  (6.2)  we  obtain  the  bound 


II Si  -  S\\  <  S(N  +  2 )u  +  2N(N  +  2)5^2||x||2u2  +  N2||*||2u2(l  +  (N  -j-  2 )u). 


Recalling  the  definition  (1.3)  of  k,  we  see  that 


<  (N  +  2 )u  4-  2W(N  +  2)«u2  4-  NVu2(1  4~  (N  4-  2 )u) 
~  Nu  4"  A/2/c2u2  4-  2 N2«u2. 


When  /V  »  1,  «  »  1,  the  term  N2k2u2  may  cause  problems  as  was  seen  in 
table  2.  Note  that  this  term  results  from  the  term  N(jj  ]T]x,£,)2  in  (6.2).  Wc  will 
now  show  that  the  computed  value  S?  is  a  good  approximation  to  this  error.  We 

have  that 
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$2  "*  ^0Ox‘  ~  *)(*  +  7.))  with  |7,|  <  (N  +  2)u  -f  0(u2). 

=  —  *)  +  (*  —  *))(!  +  7.)) 

=  Xf(S(*‘-iH—^(E**)(A,+E-w)) 

-  rO  -*>*)’  -  ^(E<*  -  **)(£  *■«.)(" + X>.) 

+  tME  (N’ + 2N  E  ->■ + (E  ■*)’) 

Note  that  S2  contains  a  term  so,  using  Bjorck's  modification  of  the 

two-pass  algorithm  we  compute  S  as 

5=  (Si—  Sy(l  + 6),  with  |«5|<u 

= (s + El*' — tf’1-  ~  (E *<&)  E<*  -  *w + *> 

-^(E(*-^(E*4W+E*) 

-^(E*#"E*+(E*)’)>+‘>- 

Bounding  these  quantities  as  before  gives  the  following  bound  for  the  relative  error: 


<  (N  -f-  2)u  -f  (N  -f-  2)('1/Vk  +  2(N-f-  2))u2 

+  (N  +  2)(3N  V  +  (G/V  -f  4)«  -f  (N  +  2 ))u3  -f-  0(u4) 


~A/u  +  4N  W  -f  3/V3kV. 


The  modification  has  thus  reduced  the  “second  order  term"  by  roughly  a  factor 
of  K. 


7.  Calculation  of  the  mean  in  double  precision. 

A  greater  accuracy  can  be  achieved  from  any  algorithm  for  computing  the 
sum  of  squares  by  simply  using  higher  precision  arithmetic.  It  is  important  to 
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note,  however,  that  a  large  increase  in  accuracy  can  often  be  achieved  by  shifting 
only  some  of  the  calculations  to  double  precision.  From  the  error  analysis  of  the 
two-pass  algorithm,  we  see  that  computing  the  sample  mean  in  double  precision 
would  replace  the  bound  |£,|  <  Nu  -f-  0(u2)  in  (6.1)  by  |£,|  <  Nit2  -j-  G(t<J).  If 
the  remainder  of  the  calculations  arc  still  computed  in  single  precision,  the  error 
bound  (6.3)  will  nonetheless  be  replaced  by  the  improved  bound 

<Nu  +  0(u3). 

The  difference  which  this  can  make  in  practice  is  evident  from  table  6  in  section 
8,  which  gives  the  results  of  some  numerical  experiments. 

The  generalized  updating  algorithm  and  the  pairwise  algorithm  arc  also 
improved  by  calculating  the  corresponding  running  sums  in  double  precision. 
Numerical  results  for  these  modified  algorithms  are  given  in  tables  7  and  8  respec¬ 
tively. 

8.  Experimental  results 

All  of  the  results  presented  in  this  section  were  computed  on  an  HIM  370/108 
computer  at  the  Stanford  Linear  Accelerator  Center.  The  data  used  was  provided 
by  a  random  number  generator  with  mean  1  and  a  variety  or  different  variances 
a2.  For  this  choice  of  the  mean,  k  \fo.  In  each  case  (he  results  have  been 
averaged  over  20  runs.  Single  precision  was  used  in  most  of  (he  u-  Is  except  in 
the  cases  where  the  mean  was  computed  in  double  precision  (table  ;  0-8).  In  single 
precision,  u  5  X  10  .  The  "correct”  answer  for  use  in  computing  (in'  error 

was  computed  in  quad  precision.  Wc  report  the  number  of  correct  digits  in  the 
calculation,  defined  as  — logoff)  where  E  is  the  relative  error. 


Si— 5 
S 
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Table  1:  Number  of  correct  digits  for  the  textbook  algorithm  on 
chosen  randomly  from  N(1.0,  a 2). 


Table  2:  Number  of  correct  digits  for  the  two-pass  algorithm  on 
chosen  randomly  from  N(1.0,  a2). 


64 

256 

1024 

2048 

1.0 

5.2 

5.1 

4.0 

4.0 

10-1 

5.4 

4.5 

4.2 

4.2 

10-2 

5.6 

4.5 

4.4 

3.7 

5.6 

4.6 

4.5 

3.6 

io~4 

5.2 

4.8 

4.4 

4.0 

10~5 

5.5 

5.3 

3.1 

3.0 

10~8 

4.5 

4.4 

2.1 

1.9 

10“7 

3.5 

3.3 

1.1 

0.9 

10-6 

2.5 

2.3 

0.1 

—0.1 

N  data  points 


N  data  points 
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Tabic  3:  Number  of  correct  digits  for  Bjorck’s  two-pass  algorithm  on  N  data  points 
chosen  randomly  from  N(1.0,  a2). 


N 


64 

256 

1024 

2048 

1.0 

5.2 

5.1 

4.0 

4.0 

10-1 

5.4 

4.5 

4.2 

4.2 

10“2 

5.6 

4.5 

4.4 

3.7 

10~3 

5.6 

4.6 

4.4 

3.6 

10— 4 

5.2 

4.8 

3.9 

3.7 

10“5 

5.2 

5.0 

3.9 

3.8 

10“6 

5.4 

5.0 

4.1 

4.0 

10-7 

5.7 

4.6 

4.2 

4.2 

10-8 

6.2 

4.6 

3.8 

3.3 

Table  4:  Number  of  correct  digits  for  the  generalized  updating  algorithm  on  1024 
data  points  chosen  randomly  from  N(1.0,  o2)  with  various  values  of  m.  (Note 
that  m  =  1  corresponds  to  algorithm  (1.5)  while  m  =  1024  is  just  the  two-pass 
algorithm). 


\  m 

1 

2 

4 

8 

16 

32 

64 

128 

256 

512 

1.0 

4.0 

4.0 

4.3 

4.6 

4.9 

5.0 

5.0 

5.1 

5.0 

4.2 

io-1 

4.2 

4.2 

4.5 

4.8 

5.0 

5.1 

5.2 

5.3 

4.5 

4.3 

1 

o 

4.5 

4.4 

4.7 

5.0 

5.1 

5.3 

5.4 

4.9 

4.5 

4.4 

10— 3 

4.1 

4.2 

4.5 

4.8 

5.0 

5.2 

5.2 

4.8 

4.6 

4.6 

10~A 

3.6 

3.7 

4.0 

4.3 

4.5 

4.8 

5.0 

4.8 

4.8 

4.6 

10~5 

3.2 

3.4 

3.8 

4.1 

4.3 

4.6 

4.8 

5.1 

5.1 

3.4 

10-6 

2.5 

3.0 

3.3 

3.7 

4.0 

4.3 

4.4 

4.6 

4.6 

2.4 

10“7 

1.4 

2.1 

2.5 

2.9 

3.6 

3.6 

3.5 

3.4 

3.3 

1.4 

10-8 

0.4 

1.0 

1.7 

2.4 

3.0 

2.8 

2.5 

2.4 

2.3 

0.4 
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Table  7:  Number  of  correct  digits  for  the  generalized  updating  algorithm  on  102-1 
data  points  chosen  randomly  from  N(i.O,  a2)  with  various  values  of  m.  In  this 
test  the  running  sums  were  computed  in  double  precision.  (Note  that  m  =  1 
corresponds  to  algorithm  (1.5)  while  m  =  1024  is  just  the  two-pass  algorithm). 


*  2  X. 

1 

2 

4 

8 

16 

32 

64 

128 

25G 

512 

1.0 

4.0 

4.0 

4.3 

4.6 

4.9 

5.0 

5.1 

5.1 

5.0 

4.2 

10-1 

4.2 

4.2 

4.5 

4.8 

5.0 

5.1 

5.2 

5.3 

4.5 

4.3 

10~2 

4.4 

4.4 

4.7 

4.9 

5.2 

5.3 

5.4 

4.9 

4.5 

4.4 

10“3 

4.4 

4.4 

4.7 

4.9 

5.2 

5.3 

5.3 

4.8 

4.0 

4.G 

10~4 

3.9 

3.9 

4.2 

4.5 

4.8 

5.0 

4.9 

4.8 

4.8 

4.8 

10-5 

3.9 

3.9 

4.2 

4.5 

4.8 

5.0 

5.0 

4.9 

4.9 

4.3 

10~6 

4.1 

4.1 

4.3 

4.G 

4.9 

5.0 

5.1 

5.1 

4.8 

4.2 

10~7 

4.2 

4.2 

4.5 

4.8 

5.0 

5.2 

5.2 

5.3 

4.5 

4.3 

10-8 

4.4 

4.4 

4.7 

5.0 

5.2 

5.3 

5.4 

4.9 

4.5 

4.4 

Table  3:  Number  of  correct  digits  for  the  pairwise  algorithm  on  N  data  points 
chosen  randomly  from  N(1.0,  a2).  In  this  test  the  running  sums  were  computed  in 
double  precision. 


N 


64 

256 

1024 

2048 

1.0 

5.9 

5.8 

5.7 

5.7 

io-1 

5.9 

5.7 

5.7 

5.7 

10-2 

6.0 

5.7 

5.7 

5.6 

10~3 

6.0 

5.7 

5.6 

5.5 

10~4 

5.9 

5.8 

5.7 

5.6 

10—5 

5.9 

5.8 

5.6 

5.6 

IO-0 

5.9 

5.8 

5.7 

5.6 

10~7 

6.0 

5.8 

5.7 

5.7 

10-8 

6.1 

5.8 

5.8 

5.6 
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t.  A  FORTRAN  implementation  o I  the  pairwiae  algorithm. 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


subroutine  updated, n, sun, s,x) 

INTEGER  N,N 
REALMS  3,SUE,X(N) 

GIVEN  THE  SUM  AND  SUM  OP  SQUARES  OP  DEVIATIONS  PROK  THE 
KEAN  FOR  A  SAMPLF  OP  M  POINTS, 


SUM 


SUM  Y  (I) 
1*1 


s  *  sun  (Y  (i) 

1*1 


2 

3UM/K) 


AND  GIVEN  K  NEW  DATA  POINTS  X(1)...X(N),  THIS  ROUTINE  PRODUCES 
THE  SUM  AND  SUM  CF  SQUARES  POS  THE  COMBINED  SAMPLE: 

N 

SU.l  :=  SUM  ♦  SUM  X(I) 

1*1 

K  2  N  2 

S  :=  SUM  (Y  (I)  -  SUM/(K*N))  ♦  SUM  (X(I)  -  SU.M/(M*N)) 

1=1  1*1 


IH2  SUM  AND  SUM  OF  SQUARES  FOR  THE  NEN  POINTS  ARE  CALCULATED 
USING  THE  PAIRWISE  ALGORTIHM.  THE  OLD  SUM  AND  SUM  OP  SQUARES 
IS  THEN  UPDATED. 

THIS  SOUTINE  HAS  LOCALLY  DIMENSIONED  ARRAYS  TERMS,  SUHA  AND 
S A  WHICH  CURRENTLY  HAVE  DIMENSION  21.  THIS  LIMITS  THE 
NUMBER  OF  POINTS  WHICH  CAN  BE  HANDLED  TO  N  <=  2**20  *  1048576. 
TO  USC  WITH  LARGER  N,  INCREASE  THESE  DIMENSIONS  TO  SOMETHING 
AT  LEAST  AS  LARGE  AS  LOG2(N)«1. 

INTEGER  TERMS  (21) , TOP, T 
R  EAL*tJ  SUM  A  (21)  ,SA  (21)  ,NEAN  ,  N  SUM  ,  NS 

TERMS  (1)  *  0 
TOP  *  2 
N 2  *  N/2 

IP  (N  .IE.  0)  GO  TC  70 
IP  (N  .GT.  1)  GO  TO  6 
NSU.M  *  X  (1) 

NS  *  0 
GO  TO  50 


on  on  on  on  non 


6  00  20  1*1, N2 

#  COMPUTE  THE  SDH  AND  SUH  OF  SQUARES  FOR  THE  NEXT  TWO 

•  DATA  POINTS  IN  X.  PUT  THESE  QUANTITIES  ON  TOP  OF 

•  THE  STACK. 

SUHA  (TOP)  *  X  (2*1-1)  ♦  X  (2*1) 

SA  (TOP)  =  (X  (2*1)  -  X  (2*1-  1)  )  **2  /  2.0 

TERNS  (TCP)  =  2 

13  IF  (IEaKS(TOP)  .NE.  TERES  (TOP- 1) )  GO  TO  20 

•  TOP  TWO  ELEMENT S  ON  STACK  CONTAIN  QUANTITIES  COMPUTED 

*  FROM  THE  SANE  NUMBER  OF  DATA  POINTS.  COMBINE  THEM: 

DP  *  TOP-1 

TERES  (TOP)  -  2  *  TERES  (TOP) 

SA  (TCP)  =  SA  (TOP)  ♦  S A  (TOP*  1)  ♦  (SUHA  (TOP)  -  SUMA  (TOP*  1)  )  **2 
X  /  TERNS  (TOP) 

SUHA  (TOP)  =  SUMA  (TOP)  ♦  SUHA(TOP*1) 

GO  TO  10 
23  TOP  =  TCP*1 

TOP  =  TOP-1 

IF  (2*N2  .EQ.  N)  GO  TO  30 

t  I  IS  ODD.  PUT  LAST  POINT  ON  STACK: 

TOP  *  TC  P* 1 
TERNS  (TOP)  =  1 
SUHA  (TOP)  *  X  (N) 

SA  (TO  P)  *  0.0 
30  T  *  TERMS  (TOP) 

NS JM  =  SUHA  (TOP) 

NS  =  SA  (TOP) 

IF  (TOP  .LT.  3)  GC  TO  50 

»  N  IS  NOT  A  POWER  OF  2,  THE  STACK  CONTAINS  MORE  THAN 

*  ONE  ELEMENT.  COMBINE  THEM: 

DO  40  J=  3, TOP 

I  *  TOP*2  -  J 

iS  *  NS  ♦  SA  (I)  ♦  T*(TBPMS  (I)  *NSUM/T  -  SUMA(I))**2  / 

X  (TERMS  (I)  *  (TERNS (I)  ♦!» ) 

NSIIM  =  NSUM  ♦  SUHA  (I) 

43  T  *  T+TERKS  (I) 


50  CONTINUE 

C  t  COMBINE  NS  AND  NSUM  WITH  S  AND  SUH  RESPECTIVELY: 

IF  (M  .2Q.  0)  GO  TO  60 

NS  =  S  «•  NS  ♦  N*  (N*SUM/E  -  NSUM)**2  /  (N*  (*«•!!)) 
NSUM  =  SUM  4  NSUM 
60  S  *  NS 

SOM  *  NSUM 
C 

70  RETURN 
END 
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